/** @file Matrix operations */
#pragma once

#include "core/geometry.hpp"

/** A 4x4 matrix */
class Matrix
{
public:

  /** The matrix data */
  float m[4][4];

  /** Check for nans matrix
   * @return if there are nans */
  inline bool hasnans()
  {
    #pragma unroll
    for(int i=0; i<4*4; i++)
      {
	if(isnan(m[int(i/4)][i%4])) return true;
      }
    return false;
  }

  /** Construct */
  Matrix() {}

  /** Construct
   * @param M matrix data */
  Matrix(float* m);

  /** Construct */
  Matrix(float c00, float c01, float c02, float c03,
	 float c10, float c11, float c12, float c13,
	 float c20, float c21, float c22, float c23,
	 float c30, float c31, float c32, float c33);

  /** Construct
   * @param o other matrix */
  Matrix(const Matrix& o);

  /** Get the transpose of the matrix
   * @return transposed matrix */
  Matrix transpose() const;

  /** Get the inverse of the matrix
   * @return inverse matrix */
  Matrix inverse() const;

  Matrix operator*(const Matrix& o) const;

  void operator*=(const Matrix& o);

  inline Vector operator*(const Vector& o) const
  {
    return Vector(o.x*m[0][0] + o.y*m[0][1] + o.z*m[0][2],
		  o.x*m[1][0] + o.y*m[1][1] + o.z*m[1][2],
		  o.x*m[2][0] + o.y*m[2][1] + o.z*m[2][2]);
  }

  inline Point operator*(const Point& o) const
  {
    return Point(o.x*m[0][0] + o.y*m[0][1] + o.z*m[0][2] + m[0][3],
		 o.x*m[1][0] + o.y*m[1][1] + o.z*m[1][2] + m[1][3],
		 o.x*m[2][0] + o.y*m[2][1] + o.z*m[2][2] + m[2][3]) / (o.x*m[3][0] + o.y*m[3][1] + o.z*m[3][2] + m[3][3]);
  }

  bool operator==(const Matrix& o) const
  {
    #pragma unroll
    for(int i=0; i<4; i++)
      {
	#pragma unroll
	for(int j=0; j<4; j++)
	  {
	    if(o.m[i][j] != m[i][j]) return false;
	  }
      }
    return true;
  }

  inline bool operator!=(const Matrix& o) const
  {
    return !(*this == o);
  }
};

inline Vector operator*(const Vector& o, const Matrix& m)
{
  return m*o;
}

inline Point operator*(const Point& o, const Matrix& m)
{
  return m*o;
}

